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ABSTRACT 

The growth of structure from scale-free initial conditions is one of the most 
important tests of cosmological simulation methods, providing a realistically complex 
problem in which numerical results can be compared to rigorous analytic scaling 
laws. Previous studies of this problem have incorporated gravitational dynamics and 
adiabatic gas dynamics, but radiative cooling, an essential element of the physics of 
galaxy formation, normally introduces a preferred timescale and therefore violates the 
conditions necessary for self-similar evolution. We show that for any specified value of 
the initial power spectrum index n [where P(k) tx k n ], there is a family of power-law 
cooling functions that preserves self-similarity by ensuring that the cooling time in an 
object of the characteristic mass M* is a fixed fraction tc of the Hubble time. We 
perform hydrodynamic numerical simulations with an Einstein-de Sitter cosmology, 
a baryon fraction of 5%, Gaussian initial conditions, two different power spectrum 
indices, and four values of tc for each index, ranging from no cooling to strong cooling. 
We restrict the numerical simulations to two dimensions in order to allow exploration 
of a wide parameter space with adequate dynamic range. In all cases, the simulations 
are remarkably successful at reproducing the analytically predicted scalings of the 
mass function of dissipated objects and the gas temperature distributions and cooled 
gas fractions in collapsed systems. While similar success with 3-D simulations must 
still be demonstrated, our results have encouraging implications for numerical studies 
of galaxy formation, indicating that simulations with resolution comparable to that in 
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many current studies can accurately follow the collapse and dissipation of baryons into 
the dense, cold systems where star formation is likely to occur. 

Subject headings: Methods: numerical — Hydrodynamics — Radiative transfer - 
Galaxies: formation — Cosmology: theory — Large scale structure of the universe 

1. Introduction 

Studies of evolution from scale-free initial conditions provide idealized but illuminating 
examples of the more general process of hierarchical structure formation. So long as the 
background cosmology, input physics, and initial conditions are scale-free, even an inherently 
complex and highly nonlinear process such as gravitationally driven, hierarchical structure 
formation must evolve self-similarly in time. Self-similar scaling offers a powerful analytic guide 
to the behavior of such complex systems, and investigations of scale-free clustering have yielded 
important insights concerning the growth of cosmological structure (e.g., Davis & Peebles 1977; 
Kaiser 1986; Efstathiou et al. 1988). The evolution of scale-free initial conditions is also one of the 
few cosmological problems in which numerical simulations can be tested against rigorous analytic 
predictions, and the study of self-similar gravitational clustering by Efstathiou et al. (1988) is one 
of the key pieces of evidence for the accuracy of cosmological N-body methods. 

In Owen et al. (1998b; hereafter Paper I), we extended this approach to adiabatic gas 
dynamics, presenting a set of smoothed particle hydrodynamics (SPH) simulations designed to 
study self-similar evolution of structure in a mixed baryon/dark matter universe. We found that 
the resulting population of collapsed structures demonstrated the expected self-similar scalings, 
so long as we were careful to account for the numerical limitations of each experiment. However, 
while the models presented in Paper I included gravitational, pressure, and shock processes, they 
neglected the effects of radiative energy loss from the gas, a critical element in the formation 
and evolution of galaxies (White & Rees 1978). In this paper we extend our earlier work to 
include radiative dissipation; for each choice of the scale-free initial power spectrum P(k) oc k n , 
we construct artificial cooling laws that maintain the scale-free nature of the physics. 

Our approach in this investigation differs from that of Paper I in a few important respects. 
First, in Paper I we considered a set of 3-D simulations, while for this study we restrict ourselves 
to 2-D simulations. The restriction to 2-D allows us to perform a larger number of high dynamic 
range experiments than would be practical in 3-D. While 3-D simulations are clearly necessary 
for realistic studies of galaxy formation, for our present purposes we wish to study the effects of 
radiative dissipation in a variety of idealized, scale-free models, and 2-D experiments provide an 
economical starting point. Our investigation represents a first attempt to examine self-similar 
evolution in models that incorporate the physical processes most essential to galaxy formation: 
gravitational collapse and merging, shock heating, and radiative cooling. We will use our results 
to guide the choice of parameters for more expensive, 3-D models in a future study. 
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Section 2 discusses the analytic scaling of characteristic group properties in 2-D models. 
In Section 3 we derive the cooling functions that preserve self-similarity, for both 2-D and 3-D 
models. Section 4 describes the simulations and their numerical limitations, and it presents 
our first important numerical result, the scaling of the mass functions of dark matter, baryon, 
and dissipated baryon groups. Section 5 presents tests of mass, temperature, density, and 
Bremsstrahlung luminosity scaling, similar to those used in Paper I. Section 6 examines the central 
issue of the paper, the scaling of the cooled baryon component in collapsed objects. In Section 
7 we summarize our conclusions and discuss their implications for numerical studies of galaxy 
formation. 



2. Expected Scalings for 2-D Perturbations 

While working in 2-D affords us a number of computational advantages, studying cosmological 
structure in 2-D introduces a number of subtle departures from the more familiar 3-D case. It is 
important to understand that these models really represent the imposition of 2-D perturbations 
in a formally 3-D universe. The correct conceptual way to view the setup for these simulations is 
that we are taking an initially homogeneous, 3-D universe containing baryons and dark matter and 
initializing only k x and k y density perturbations in Fourier space. These initial conditions imply 
that structure evolves in the x and y directions alone - matter continues to simply expand with 
the Hubble flow in the z direction. Therefore, the most collapsed structures that form represent 
infinitely long filaments, rather than collapsed groups or galaxies. When we measure the "mass" 
of one of these 2-D objects, we are in fact measuring a mass per unit length, or linear density 
9 oc m/t. Similarly, the mass of each particle in the simulations is actually a mass per unit length, 
and the particles interact gravitationally as though they were infinite, thin rods. 

To illustrate some of the more important distinctions from the general 3-D case, consider the 
collapse and formation of an isolated filament. If the total mass per unit length of this object is 
6, and we have another infinite, thin rod of mass per unit length at a distance r, then their 
mutual gravitational force is F = —G9@/r oc 1/r, and the gravitational potential of the first rod 
is <J>(r) = —GOlnr. The circular velocity, v 2 c = rd^/dr = GO, is independent of distance. In order 
to define a temperature that corresponds to this circular velocity (see Thoul & Weinberg 1996 
for the 3-D analogy), consider an isothermal atmosphere of gas about the filament in hydrostatic 
equilibrium. If the gas follows a density profile p(r) = p\(r jr\)~ a and the pressure is given by 
p(r) = fcsT(//m p ) _1 /9(r), then the requirement of hydrostatic equilibrium becomes 

GO ^ . , . dp n 

— 2irr ar pir) = — — 2irr ar, (1) 
r dr 



which after some manipulation yields 



T = ^G9 = ^Evl (2) 
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Inspection reveals that typically our 2-D objects obey p(r) oc r 2 (where p is the mass per unit 
volume), so we recover the 3-D result, T = pm p (2kB)~ 1 v 2 . 

If the input physics, background cosmology, and initial conditions are scale-free, then at a 
given time the only characteristic physical scale is set by the amplitude of density fluctuations. 
This condition implies scaling laws for the time evolution of radii, masses, temperatures, and 
Bremsstrahlung luminosities of collapsed objects. We derived these scaling laws in Paper I, 
following the reasoning of Kaiser (1986), but here we must reformulate them for 2-D perturbations. 
The variance of (linear theory) mass fluctuations smoothed with a window of comoving scale R c is 

a 2 M (R c ,a) = / 2TrkdkP(k,a)W(kR c ) oc a 2 (i? c )- (n+2) , (3) 
Jo 

where P(k,a) oc a 2 k n is the linear fluctuation power spectrum. We can define a characteristic 
comoving length scale R1 to be the scale on which the linear fluctuation variance has some 
specified value. Equation (Q) implies that R% oc a 2 ^ 2+n \ and the corresponding radius in physical 
units is 

R* ocaRt oca( 4+n )/( 2+ "). (4) 
Since the underlying physics is 3-D, the proper background density scales as a -3 , and we have 

oc p oc a~ 3 . (5) 

The characteristic linear density 9* (analogous to the characteristic mass M* in 3-D) scales as 

9<; o, pt(R c *) 2 oc a^ 2+n > ^oca-X c oca( 2 - n )/P+"). (6) 

The difference between the physical 6* and the comoving 9% is the a -1 factor, which reflects the 
drop of linear mass densities caused by expansion along the 3rd axis. Based on equation @, we 
know that the hydrostatic equilibrium temperature scales directly with 9, so 

T* oc0* oc a {2 -' n)/{2+n \ (7) 

Finally, the Bremsstrahlung luminosity per unit length £ oc L/£ scales as 

£oc^T, 1/2 oca-( 6+9 ")/( 4+2 "). (8) 

Note that we are defining the luminosity to be proportional to the square-root of the temperature. 
Even though we will be assuming a different temperature dependence for radiative cooling in the 
following section, whenever we refer to the "luminosity" we will still use C oc T 1 / 2 , since this 
corresponds to the physical process of Bremsstrahlung emission in the real universe. 

The *-ed quantities in the proportionality relations could refer to any specific choice 

of the linear fluctuation variance. In §4 we will give a precise definition of 9* in terms of the 
characteristic mass appearing in the Press-Schechter (1974) mass function formula (eq. |l^] below). 
We will also refer frequently to the mass scale 9 n \ on which the fluctuation variance is unity, 
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indicating that rms fluctuations on this scale are entering the non-linear regime. Specifically, this 
mass scale is defined by the conditions 

e nl (a)=7rR 2 nl p(a), (9) 

where 

a 2 (R nh a) = l. (10) 



3. Radiative Cooling Laws 

A realistic cooling function like the one arising in a primordial H/He plasma would impose a 
dimensional physical scale into our calculations, violating the conditions that lead to self-similar 
time evolution. The clearest demonstration of this problem is that the cooling time in an object 
of characteristic mass M* would not be a constant fraction of the Hubble time, since the standard 
cooling function has a complicated temperature dependence and the characteristic temperature 
T* and density p* evolve in time. For example, at some epochs most of the baryonic mass in an 
M* object might be able to cool in a Hubble time, while at other epochs only a small fraction 
would be able to cool. We therefore cannot use the standard physical cooling law and maintain 
the scale-free nature of the problem. 

Nonetheless, for a given power-law initial power spectrum, P(k) oc k n , it is possible to 
determine a power-law cooling relation A(T)/nn 2 oc T 13 such that the cooling time in an M* 
object is a fixed fraction of the Hubble time. So long as the cooling law meets this requirement, 
the physical system remains scale-free, and the rigorous prediction of self-similar evolution for the 
system holds. Under this condition, at a fixed expansion factor the ratio of the cooling time to the 
Hubble time may vary with mass, but the ratio depends only on M/M* and thus scales properly 
with time. The dependence of the cooling exponent (5 on the power spectrum index n is different 
for 2-D and 3-D perturbations, so we will give the results for both cases here. 

The radiative contribution to the time rate of change of the specific thermal energy u is 



Du 
~Dt 



-42. (id 

R P 



where we use the subscript R to denote the fact that this is only the radiative piece of Du/ Dt. 
We parameterize the radiative cooling function A(T) as a power-law, 

= [pm p (l + n Hc /n H )] 2 ^ = A T? = A x u? . (12) 

«h p 

In these terms the cooling time tc can be expressed as 



tc = « ( ^ I J _1 = « (^P) 1 = u( P A lU ^ = A^p-W-P. 



(13) 
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AT'Po 1 ^ = i c H, 1 ( ± ) . (14) 



The requirement to maintain self-similarity is that the cooling time t c for a characteristic 9* (or 
M*) object must be a fixed fraction £c of the Hubble time i# = i?5~ 1 (a/ao) 3//2 . Note that in 
general we expect the cooling time to be a function of mass — e.g., tc(0.5 0*) ^ to(2 0*) — but 
we are forcing the cooling time for any given multiple of a characteristic mass to be fixed at all 
times, so that tc(2 di) = tc(2 6*, 02). 

For 2-D perturbations we know that the specific thermal energy and the density scale as 
u* oc a( 2 ~ n )' / ( 2+n ) and oc a~ 3 , so fixing = ictu yields, after some manipulation, 

-3/2-(l-/3)(2-n)/(2+n) 

Our requirement that t* c remain the same fraction of the Hubble time at all times implies that the 
a dependence on the right-hand side of equation (O) must vanish, which finally gives us 

02-D = 1^ + 1- (15) 
2 2 — n 

The 3-D derivation is very similar, except that for 3-D perturbations the specific thermal energy 
scales as u* oc ( 3+ri ) . Following the same arguments as above yields 

Now that we know the appropriate index for the cooling power-law, we need only specify the 
normalization. According to equation ([14]), the normalization A% is proportional to 

M oc ^ oc -3^, (17) 

where po and To are chosen as characteristic of a fiducial group at expansion a = a,Q. We will adopt 
the convention that the characteristic density is po ~ lOOOpbary) a s this is typical of the objects 
we select. A well-defined choice for the characteristic temperature is the hydrostatic equilibrium 
support temperature for an object with the nonlinear mass, 

T n i = T hc (0n\) = pm p {2k B )~ l Ge^. (18) 

Then by choosing a cooling time fraction tc = tc/tH, we completely specify the normalization of 
the cooling law. For tc <C 1, most of the baryon mass in 9 n \ objects should be able to cool, while 
for t c ^> 1 most of the baryons in a 6 n \ object will be unable to cool before the it merges with 
another of comparable mass, shock heating the gas to a higher temperature. 



4. The Simulations 



We examine a number of different physical scenarios, based upon two distinct sets of initial 
density perturbations and various choices for the amplitude of the associated cooling law. As in 
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Paper I, the background cosmology is a flat, Einstein-de Sitter universe with f^bary = 0.05 and 
^dm = 0.95. We consider two sets of initial density perturbations, each Gaussian distributed, with 
power-law power spectra of n = and n = +1 [where P(k) oc k n ]. Dynamically these correspond 
to 3-D perturbation spectra of indices ri3_D = — 1 and 0, respectively, for quantities represented by 
integrals over the power spectrum (such as the rate at which the scale of nonlinearity grows with 
time). 

The initial amplitude of the density fluctuations is chosen so that the linearly extrapolated 
mass fluctuations at the end of the simulation have an rms amplitude of unity for a top-hat filter 
of radius i?th ~ 0.1Lb ox for n = and Rth ~ 0.3Lb ox for n = +1, where Lbo X is the simulation 
box size. Throughout this paper we parameterize the evolution in terms of the expansion factor 
a, defined so that the final output expansion in each scenario is at = 1. In principle, because 
of the scale-free nature of these experiments, there are many possible choices for identifying the 
simulation box size with a physical scale in the real universe. For instance, if we identify the 
a,f = 1 output with the observed universe at z = and adopt the normalization condition erg = 0.5 
(appropriate for the rms mass fluctuation in spheres of radius 8 h^ 1 Mpc with Q = 1; see White, 
Efstathiou, & Frenk 1993), then the implied comoving box size is L^ ox ss 40 h" 1 Mpc for the 
n = models, and -Lbox ~ 17 hr x Mpc for n = +1. Earlier output times can be identified with 
higher redshifts, z = a -1 — 1. Alternatively, one can identify any expansion a with the z = 
universe, in which case the implied box size (for <7g = 0.5) is Lb ox = 40 a -1 h^ 1 Mpc for n = 0, 
and Lbox = 17 a~ 2 / 3 h~ 1 Mpc for n = +1, decreasing steadily as the nonlinear scale becomes a 
larger fraction of the box size. 

For both the n = and n = +1 initial conditions, we perform simulations with several 
amplitudes of the cooling law, parameterized in terms of the dimensionless cooling time 
tc = tc/tn as tc = oo (i.e., no cooling), 200, 1, and 0.1. All simulations are performed using 
•^Vbary = Ad m = 128 2 computational nodes, except for one case where we repeat the n = 0, tc = 1 
model with Ab ar y = ^dm = 256 2 nodes in order to check for convergence. The gravitational 
interactions are evaluated on a 512 2 Particle-Mesh (PM) grid, and the hydro dynamical interactions 
are modeled using Adaptive Smoothed Particle Hydrodynamics (ASPH). For a complete 
description of the code and techniques used, see Owen et al. (1998a). In all we present nine 
simulations of eight distinct physical models (2 sets of initial density perturbations x 4 possible 
cooling law amplitudes). 

We examine each simulation at expansions spaced logarithmically, roughly with an interval 
of A log a = 0.2. We are primarily interested in the properties of the collapsed objects that form 
in each model, particularly in those that undergo significant radiative dissipation. Therefore, 
as in Paper I, at each expansion we identify groups of particles using the "friends-of-friends" 
algorithm (see, e.g., Barnes et al. 1985; the specific implementation used here can be found 
at http: / / www-hpcc.astro.washington.edu / tools /FOF/| ) . We compute global average 



properties for each group (such as the dark matter and baryon mass, the mass and emission 
weighted temperatures, luminosity, and so forth) and check for self-similar evolution in the 
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distributions of properties of these groups. Although there are more sophisticated group finding 
algorithms available, friends-of-friends provides a simple, flexible, and unambiguous definition of a 
group using an algorithm that maintains the conditions necessary for self-similar scaling (because 
it does not introduce a fixed physical scale). In analyzing these experiments we use linking 
parameters of £ = 0.2 Ax p (the same as used in Paper I) and £ = 0.05 Ax p . Groups selected with 
£ = 0.2 Ax p consist of both hot and dissipated gas, whereas £ = 0.05 Ax p tends to select objects 
that have undergone significant radiative dissipation. There are essentially no dark matter objects 
found with £ = 0.05 Ax p , and very little baryonic mass collapses to this density in the models 
without cooling. 
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Fig. 1. — Comoving 9 resolution limits as a function of expansion for the (a) n = and (b) n = +1 
simulations. Although oc m/I is a linear mass density, in this 2-D context it plays the same 
role as the mass in 3-D. The horizontal dotted lines show the low 6 limit, assumed to be 28 times 
the particle "mass" for the iV bary = A^m = 128 2 simulations (actually the linear mass density per 
particle: 6i = L^ ox p/N). The dashed lines show various multiples of the nonlinear mass scale 6 n \, 
while the solid lines show the PS prediction such that statistically we would expect to find 1 (#i), 
10 (6*io), and 100 (0ioo) objects of that mass in the simulation volume. All 0's are scaled such that 
the total 9 in the volume is 1. 

In order to understand the regimes we can probe, we must first identify the mass resolution 
limits of our experiments. Since ASPH is a Lagrangian technique, the lower limit on the 
hydrodynamic interactions is best expressed as a mass limit, set by a multiple of the particle mass. 
In each experiment the ASPH smoothing tensors are evolved so that each node sees significant 
contributions from roughly 28 of its neighbors, which provides a reasonable lower limit on the 
ASPH mass resolution. The gravitational force resolution is effectively determined by a multiple 
of the gravitational softening length, in this case the size of a PM grid cell, so unfortunately this 
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is not Lagrangian in the same way as the hydro dynamic resolution. However, the PM grid has 
a linear resolution of 512 elements across the system, and on the scale of the objects we will be 
examining we may consider the gravitational resolution to be unrestrictive. As a check of this 
assumption, we have experimentally verified that the results are insensitive to changes in the 
gravitational PM resolution by factors of two (a factor of 4 in the total number of PM cells). 

The upper limit on the mass range we are sensitive to is set by our box size. The larger a group 
is, the more statistically rare it is, and the less likely we are to find examples of such structures in 
any finite, randomly realized volume. We can therefore only expect to find objects up to a certain 
size at any given output time within our finite simulation volume. The Press-Schechter mass 
function (Press & Schechter 1974; hereafter PS) provides a rough estimate of this limit. However, 
there is a subtle distinction between using the 2-D version of PS theory (eq. |l^] below) and using 
the familiar, 3-D PS formalism. In 3-D, friends-of-friends with a linking parameter i = 0.2Ax p 
selects objects within an overdensity contour of roughly 5 = 5p/p ~ 250. The well-known solution 
for the collapse of an isolated 3-D perturbation (the "top-hat" collapse model) tells us that an 
object with an actual overdensity of 5 ~ 250 corresponds to a region in the initial conditions with a 
linearly extrapolated mass fluctuation of <5* = 1.68. We can therefore use <5* = 1.68 in the PS mass 
distribution to predict the properties of the population of objects identified by friends-of-friends. 
In 2-D, however, the familiar top-hat solution does not apply, so we do not have an a priori 
prediction for the appropriate value of o* to plug into the PS mass function. Instead, we measure 
the differential group mass distribution f(9) in the simulations and fit the PS model by varying 
the value of 5*, the results of which can be seen in Figure |2| below. We find that values of 5* ~ 0.7 
for n = and <5* « 3 for n = +1, are appropriate. 

With this knowledge, we can plot a set of lower and upper comoving 9 limits for each set of 
initial conditions, as shown in Figure [I]. This figure also shows the evolution of various multiples 
of the comoving nonlinear "mass" 9^, defined earlier in equation (||). Based on these results, for 
n = we can expect most of the mass range 9 C G [1 9 n \, 10 9 n \] to be accessible over a range of 
expansions logo/a/ G [—0.7, —0.1]. Similarly, for n = +1 we find 9 C G [0.1 9 n \, 1 9 n \] should be 
accessible over log a/a/ G [—1,0]. 

In Figure |2| we plot the mass distribution f(9) of objects identified by friends-of-friends over 
a range of expansion factors for the simulations without cooling (tc = oo) and for the cooling 
models with tc = 1 and tc = 0.1. Note that f{9) is the amount of mass contained in groups 
in each mass range, not the number of objects, and that we plot the baryon and dark matter 
groupings separately. At each expansion factor f(9) is calculated and shifted (assuming that the 
self-similar relation 9 oc a'- 2_n ^'- 2+n - ) holds) to a = aj for comparison, giving the heavy solid lines. 
If the simulations perfectly obeyed the analytically predicted self-similar scalings, then the plotted 
mass functions would connect together into a continuous relation. There are no free parameters 
and no approximations in this scaling, and, as with the adiabatic results in Paper I, the measured 
distributions f(9) at different expansions link up nicely. Therefore, it appears that the group mass 
distribution is following the expected self-similar behavior well. The only notable exception is the 
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Fig. 2. — Differential mass distribution function f{6) for the (a) n = and (b) n = +1 simulations, 
scaled to a = at. In each figure the left column shows the dark matter groups linked with 
£ = 0.2 Ax p , the middle column the baryon groups linked with £ = 0.2 Ax p , and the right 
column the baryon groups linked with £ = 0.05 Ax p . The heavy solid lines show the measured 
f(9) at various expansions, shifted (assuming self-similarity) to a = at. The thin solid lines show 
the PS mass function (eq. jl9| ) assuming (a) 5* = 0.7 and (b) 5* = 3, fitted as the appropriate 
comparison for objects identified with linking parameter £ = 0.2 Ax p for the n = and n = +1 
models, respectively. Similarly, the thin dotted lines in the £ = 0.05 Ax p panels show the PS mass 
function using (a) 5* = 1 (n = 0) and (b) (5* = 4 (n = +1). 
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£ = 0.05 Ax p groups in the models without cooling (the upper right panels). This failure is to 
be expected, though, since very little mass in the non-dissipative models collapses to such high 
density, and therefore the statistical sampling in this regime is poor. 

The thin lines in Figure [2] show the 2-D PS prediction, given by the relation 

where f(9) d6 represents the mass fraction of groups in the mass range [9,9 + d0]. Here 9* is 
defined similarly to 9 n \ in equations @ and (0), except that the variance in equation (|i~0| ) is set 
to $1 instead of to unity. Altering 5* slides the resulting PS curve horizontally in this Figure but 
does not change its shape. Since, as mentioned above, we do not have the 3-D spherical collapse 
model as a guide, we fit the PS prediction to the measured f{9) curves with 5* as a free parameter, 
obtaining 5* = 0.7 for n = and <5* = 3 for n = +1. Equation (|3|) shows that the relation between 
9* and our other fiducial "mass" 9 n \ is 

9* = 9 nl 5i /{2+n) . (20) 

For the £ = 0.2 Ax p objects, PS curves with the same value of 5* fit the dark matter and 
baryonic objects equally well. This simultaneous fit implies that the baryon-to-dark matter ratio 
in these structures is near the universal average, in contrast to the results of Paper I, where we 
found that the baryons in groups tended to be underrepresented relative to the dark matter. The 
difference might reflect either the 2-D geometry or the somewhat higher resolution of our present 
experiments. There is also a hint in Figure ^ (and in Figure |3| below) that the baryon fractions are 
slightly higher in the models with radiative cooling. 

The thin dotted lines in the right columns of Figure ^ show the PS prediction using 5* = 1 
for n = and 5* = 4 for n = +1. These dotted lines match the distribution of mass in dissipated 
gas objects (as selected by t = 0.05 Ax p ) reasonably well, despite the fact that there is no rigorous 
justification for using PS theory to make predictions about such dissipative structures — PS 
theory is built on an approximate description of purely gravitational dynamics. Since there is a 
free parameter in our PS fits and it is chosen separately for £ = 0.05Ax p and i = 0.2Ax p , it is not 
clear how much significance should be attached to this match. It does appear, however, that the 
relation between "dissipated mass" (identified by £ = 0.05Ax p ) and "virialized mass" (identified 
by £ = 0.2Ax p ) is simple enough in these scale-free models that PS theory gives a reasonable guide 
to the mass function of the cooled objects. 



5. Testing for Self-Similar Evolution 

Figure ^ provides the first evidence that our simulations are following the expected self-similar 
scaling laws, at least with regard to the distributions of masses of collapsed objects. We now 
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Fig. 3. — Evolution of the fiducial comoving linear mass density 0j % for the (a) n = and (b) 
n = +1 models, where #£ 0% is defined so that at each output time 70% of the mass in the simulation 
is contained in groups with 9 C < Oj^. The points show the measured simulation results at each 
output time, while the lines show the expected scalings (9 C oc a 2 for n = and 9 C oc a 4 / 3 for 
n = +1) normalized to the data. 



turn to tests that directly examine the scaling of characteristic group masses, temperatures, 
Bremsstrahlung luminosities, and densities. 

Figure |3] shows the evolution of a fiducial group mass #£ 0% , defined so that 70% of the 
total mass is contained in groups with 9 C < 0j O y o - The "c" superscript denotes comoving linear 
mass densities; we remove the a -1 scaling so that the total "mass" 8% ot of each species in the 
simulation box is independent of time. Isolated particles are counted as "groups" with mass 
^particle f° r this purpose. We plot the four cooling models (tc = oo, 200, 1, and 0.1) for each 
initial spectral index n, simulated at our standard N = 2 x 128 2 particle resolution, and also a 
high-resolution (N = 2 x 256 2 particles) simulation of the n = 0, tc = 1 model. The lines show 
the self-similar scaling solutions (a) 9j % oc o 2 for n = 0, and (b) 6^ 0% a a 4 / 3 for n = +1 (eq. ||), 
normalized to the average amplitude of each experiment. This normalization is determined by first 
scaling each measurement to a fiducial time (assuming self-similarity), then taking the average of 
these scaled measurements. We plot both the dark matter and the baryon groups selected with 
£ = 0.2 Ax p and the baryon groups only for I = 0.05 Ax p . There are virtually no dark matter 
groups identifiable with such a small linking parameter because dissipation is required to achieve 
the necessary overdensity, at least given our finite mass resolution. 

While this Figure (and the succeeding similarly defined figures) are comparable to the 3-D 
versions in Paper I, there is one important computational distinction in how the mass fractions 
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are calculated here. In the 3-D simulations of Paper I, we typically have between several hundred 
and a few thousand identified objects, enough to define a smooth and well-behaved cumulative 
mass function. However, in the 2-D simulations we typically have an order of magnitude fewer 
objects (because we have 16 times fewer particles). These small numbers make the raw 2-D 
cumulative mass functions much jumpier than those in the 3-D simulations, which can make the 
measured #y % noisy if we only use information from those points in the cumulative mass function 
that bracket the 70% cutoff. In order to circumvent this problem, we first fit a polynomial to the 
cumulative mass function using general linear least-squares, then interpolate to find the desired 
mass fraction cutoff in this smooth polynomial fit. This method allows us to use information from 
the entire cumulative mass function in determining the 70% cutoff, and we thereby suppress the 
noise caused by the relatively small number of collapsed objects. 

Clearly all of the models demonstrate good scaling of 0j % over a range log a/a/ G [—0.8, 0] for 
n = and log a/a/ S [—1.4,0] for n = +1, similar to the ranges we expect based on the resolution 
limits plotted in Figure [l]. The scaling holds equally well for objects identified with £ = 0.2 Ax p 
(dark matter and baryons) and for dissipated baryon gas objects selected with £ = 0.05 Ax p . We 
can see some evidence that the scalings are poorer at early times, most likely due to the relative 
lack of well resolved objects at those times. The standard and high resolution versions of the 
n = 0, tc = 1 model agree almost perfectly, except that the high resolution simulation is able to 
probe back to earlier times. In both the n = and n = +1 cases, the baryon objects selected with 
£ = 0.2 Ax p in the models with cooling are slightly more massive than those in the non-cooling 
models, suggesting that dissipation plays some role in concentrating the baryons even at this 
moderate overdensity. 

In Figure |] we plot the evolution of a fiducial emission weighted group temperature T m %, 
defined in a similar fashion to the the fiducial mass used in Figure ||. The emission weighted 
temperature is defined assuming Bremsstrahlung radiation, for which the bolometric, volume 
emissivity goes as e oc p^T 1 ' 2 . This definition is somewhat inconsistent for the cases with radiative 
cooling, since the radiative energy loss from the gas follows our artificial cooling laws rather than 
a T 1//2 law. We make this choice because the Bremsstrahlung weighted temperature is closer to 
an observationally relevant quantity and because it makes it easier to compare the simulations to 
each other and to the 3-D simulations presented in Paper I. The emission weighted temperature 
associated with group i is 

feTdA _ fp*T^dA _ E 3 e jP3 Tf 2 
1 JedA Jp^/ 2 dA YjWj 

represented as a sum over the particles j that are members of group i. We sort the individual group 
temperatures in ascending order, and define the fiducial temperature T 70 % such that 70% of the 
mass is contained in objects with T < T-?q% (using the same polynomial interpolation procedure 
applied to the fiducial mass calculation). Consistent with the convention adopted for the mass 
scaling test shown in Figure ||, we count the mass contained in unresolved groups (or individual 
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Fig. 4. — Evolution of the emission weighted group temperature T 70 % for the n = (left column) 
and n = +1 (right column) models. T 70 % is defined so that at each expansion factor 70% of the 
mass in the simulation is contained in groups with T < T-jq%. The lines show the expected scalings 
normalized to the data: T oc a 1 for n = and T cx a 1 / 3 for n = +1. Point and line types as in 
Figure |3[ open squares/solid lines represent the models without cooling, stellated triangles/dotted 
lines models with cooling time tc = 200, open triangles/short dashes tc = 1, circles/long dashes 
tc = 0.1, and filled triangles/dot-dash lines the N = 2 x 256 2 , tc = 1 model. 
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particles) as having temperatures less than the lowest measured group temperature. The lines 
show the self-similar solution normalized to each experiment. Temperatures are scaled to the 
hydrostatic equilibrium support temperature (eq. 0) of a 6* object at the final expansion factor 
for each experiment: T*(aj) = G[im p (2kB)~ 1 0*(af), where 0*(df) is the mass scale associated 
with the density contrast 5* that fits the £ = 0.2 Ax p mass function, as described in §|4|. Note that 
if we were to scale instead to the characteristic temperature at each expansion factor T*(a), then 
the predicted evolution tracks would be horizontal lines. 

In Figure f| we see that the fiducial emission weighted temperature scales well in nearly all 
of the simulations, both for £ = 0.2 Ax p and for £ = 0.05 Ax p . The logarithmic slope of the 
scaling relation is three times higher for the n = models than for n = +1 (T oc a 1 vs. T oc a 1 / 3 ), 
and the simulations clearly capture this distinction (note, however, that we adapt the range of 
the horizontal axis to the dynamic range of the experiments). The standard and high resolution 
simulations of the n = 0, tc = 1 model yield nearly identical results. 

Figure |] clearly shows the effects of radiative cooling on the typical gas temperatures: the 
collapsed gas in the tc = 0.1, n = model, for example, is 1.5 dex cooler than in the case 
without cooling. Additionally, for the models with cooling, the dense, dissipated gas selected 
with £ = 0.05 Ax p is typically 0.5 dex cooler than the I = 0.2 Ax p gas. As noted in the mass 
scaling tests, we again see poorer scaling at early times, particularly noticeable in the n = 0, 
£ = 0.2 Ax p objects, where we only find reasonable scaling for log a/a/ € [—0.6, 0]. The dissipated, 
£ = 0.05 Ax p objects do not even form until the larger £ = 0.2 Ax p structures begin to scale 
well, but once this happens the temperatures of the dissipated groups always scale well. For the 
n = +1 models, the temperatures of the dissipated structures in the lower panel always scale more 
effectively than those of the larger i = 0.2 Ax p groupings. This difference suggests that the cooled 
gas cores of the structures in the n = +1 models scale correctly, and that it is the hot, diffuse gas 
in the outer regions that does not. This result seems contrary to what we find in the 3-D models 
without cooling in Paper I, and it could be due to an increased importance of hot gas that is not 
physically bound being erroneously included in collapsed structures by friends-of-friends with the 
larger linking parameter. We will consider this question further in §|(| 

Figure || tests the scaling of a fiducial group luminosity £70%! defined in much the same 

manner as the fiducial group temperature in Figure ||. For each group the total luminosity is 

defined as £ = / p 2 T x / 2 dV = J2j QjPjTj' ' ■ As with the temperature, we sort the groups in order 

of increasing luminosity and find the luminosity £70% such that 70% of the total baryon mass is 

contained in groups with £ < £70%. The luminosities are scaled to the characteristic luminosity 

1/2 

at the final expansion factor for each model, £*(a/) = 0*(a/) p*(af) T* (o/)) where we again 
use quantities appropriate for the £ = 0.2 Ax p groups. The lines show the self-similar predictions 
£ oc a~ 3 / 2 for n = and £ oc a~ 5 / 2 for n = +1 (eq. ||). The luminosity does not depend strongly 
on the adopted cooling law because the reduction in temperature in models with stronger cooling 
is compensated by the increase in gas density. 
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Fig. 5. — Evolution of the group luminosity C 70 % for the (a) n = and (b) n = +1 models, such 
that at each expansion 70% of the mass in the simulation is contained in groups with £ < jC 70 y . 
Note that the "luminosity", C, is in fact a luminosity per unit length. The lines show the expected 
scalings normalized to the data: C oc a -3 / 2 for n = and C oc a -5 / 2 for n = +1. 



Relative to the poor scaling of group luminosities found for 3-D objects in Paper I, the scaling 
behavior in Figure |5] is remarkably good. The £ = 0.2Ax p groups follow the analytic scaling laws 
for all combinations of n and tc- The dissipated, t = 0.05 groups show some tendency to have 
lower luminosities at earlier times (when the groups are less well resolved), but even in these cases 
the luminosity scaling is not bad. Furthermore, the high and low resolution versions of the n = 0, 
tc = 1 model concur almost perfectly, in marked contrast to the strong resolution dependence 
found for group luminosities in the 3-D non-cooling models of Paper I. We will return to consider 
the difference between our present and previous results, following a discussion of the gas density 
scaling tests. 

We cannot test the gas density scaling in precisely the same manner as the previous quantities 
because groups are effectively selected by their overdensities, and we might therefore expect their 
average densities to scale by construction. We approach this problem as we did in Paper I, by 
first calculating a percentile density p x % for each group, such that x% of the gas particles in the 
group are at densities pj < p x % . We then take all of the groups in a given mass range and find the 
average value of p x %, where the average is weighted by the group mass, with the constraint that 
a group must contain at least 28 particles to contribute. In Figure |(| we plot average densities 
defined in this way for p\o% , />5o% , and P9o% > progressively probing from the outskirts of each group 
inward. The left columns are calculated for objects selected with the larger linking parameter 
t = 0.2 Ax p (including both the hot and cold components of collapsed structures), while the right 
columns use the smaller linking parameter i = 0.05 Ax p (effectively restricting the selection to 
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Fig. 6. — Evolution of the group density p x % for the (a) n = and (b) n = +1 models. p x % is 
defined as the average of p x% for groups in a given range of 9, where p x% for each group i is the 
density such that x% of the mass in the group is at densities p < p x% . In each figure, the left 
column shows the results for groups selected with linking parameter I = 0.2 Ax p and the right 
column for groups selected with i = 0.05 Ax p . The rows show the averages for pio%, /?5o%> and 
/°90%; progressively sampling from the outskirts of each group inward toward the core. For the 
n = models in part (a) we average over objects in the mass range 0.5 n \ < 9 < 2.5 8 n \, while 
for the n = +1 models in part (b) we average over groups in the range 0.2 9 n \ < 9 < 1.0 9 n \. The 
lines show the expected scaling normalized to the data: p oc a" 3 . Point and line types as defined in 
Figure |3[ open squares/solid lines represent the models without cooling, stellated triangles/dotted 
lines models with cooling time tc = 200, open triangles/short dashes tc = 1, circles/long dashes 
ic = 0.1, and filled triangles/dot-dash lines the N = 2 x 256 2 , to = 1 model. 
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only the cooled gas cores in the simulations with radiative cooling). 

Considering first the t = 0.2 Ax p groups in the left columns, the effects of cooling are quite 
evident, as models with stronger cooling show markedly higher densities. We can also see that 
the behavior of the density noted in Paper I is qualitatively reproduced here (most notably in 
the n = models): the lower density cuts (probing the outer regions of the objects) scale most 
effectively, but as we consider progressively higher density cuts (and therefore move inward to 
smaller radii) the gas density progressively scales more poorly. For n = +1 the density scaling 
is better, and the differences in the scaling behavior of p\o% and p$o% are smaller, though they 
follow the same trend. The density evolution for the dissipated baryon structures selected with 
£ = 0.05 Ax p (the right columns in Figure ||) is similar to pgo% in the I = 0.2 Ax p objects and 
demonstrates relatively little change between Piq%, P^q%, and p$Q%. Presumably this similarity 
arises because the dissipated groupings selected with I = 0.05 Ax p are formed at the cores of the 
larger, less dense structures identified by t = 0.2 Ax p , and therefore pgo% in £ = 0.2 Ax p probes 
the same material selected with the smaller linking parameter. 

Comparing these 2-D results to the 3-D results of Paper I, we find that the core densities 
G°90%) scale more accurately in 2-D but that the halo densities (p\o%) scale somewhat less 
accurately. This difference most likely reflects the interplay between the dimensionality of the 
models and the Lagrangian resolution of SPH, in which quantities are averaged over fixed numbers 
of neighbors. In both 2-D and 3-D simulations we find that the average radial density profiles 
approximately follow p(r) cx r~ 2 over the scales we resolve. In 2-D, the number of particles 
in a radial shell goes as N2-d{v) cx p{r) dA = 2irrp(r) dr oc r , while in 3-D this becomes 
N(r) oc p(r) dV = Ai[r 2 p(r) dr oc r . Thus, for 2-D objects we tend to have relatively poorer 
resolution for large r (the region probed by P\o%), and better resolution for small r (pgo%)- 

The better scaling of core gas densities in these 2-D experiments partly accounts for the 
improved scaling of the luminosities relative to their 3-D counterparts in Paper I. However, the 
core gas densities do not scale perfectly, and this change alone seems unlikely to account for so 
much improvement in the luminosity scaling. The change from traditional SPH to ASPH also 
helps — we have repeated some of the 2-D tests with SPH and find that the luminosity scaling is 
somewhat worse than with ASPH but is still markedly better than that found in Paper I. The 
primary cause of the improved scaling, however, is that the luminosities of the 2-D groups are less 
core dominated and are therefore less sensitive to poorly resolved gas. The 2-D groups have nearly 
isothermal temperature profiles, while the 3-D temperature profiles fall with radius (roughly 
following T oc r -0 ' 6 ). Since the 2-D and 3-D density profiles are nearly the same (p oc r~ 2 ), the 
flatter temperature profile in 2-D leads to a less sharply peaked luminosity profile. 
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6. Scaling in the Cooled Baryon Component 

Of all the questions we can ask of our simulations, the one that goes to the heart of the 
galaxy formation issue is this: do they correctly follow the collapse and dissipative settling of 
baryons into dense, radiatively cooled objects? The scaling of the mass function (Figure |2|) and 
other properties of high density (£ = 0.05Ax p ) groups already suggests that the answer to this 
question is at least a qualified yes. In this Section, we focus more directly on the cooled baryon 
component of virialized groups. 

We begin by running friends-of-friends with linking parameter £ = 0.2 Ax p on all particles in a 
given experiment, yielding structures comprised of both baryons and dark matter (while up to now 
we have considered baryon and dark matter groups separately). Equation (||) gives a characteristic 
support temperature for an object based purely on its mass and the assumption of hydrostatic 
equilibrium: The = {jLm p (2kB)~ 1 G9. In Figure [5], we plot average baryon temperature distribution 
functions for objects that fall within a characteristic mass range, where the temperatures are 
scaled to this characteristic temperature for each structure. The average profiles are determined 
by first computing individual distributions for the baryon mass /j(T/Th e ) for each group i, then 
averaging these individual profiles in a restricted range of 8, chosen so that we will have a number 
of resolved objects in that mass range at the times we consider. 

Each panel of Figure [7| shows average temperature distributions calculated at two different 
expansion factors (a/aj = 0.4 and 1). Since the temperatures are scaled by the group hydrostatic 
equilibrium temperature and the mass range of groups contributing to the average distribution 
is scaled to 9 n \, the two curves should lie on top of one another if the simulations perfectly obey 
the expected self-similar scalings. The agreement between the dotted and solid curves in Figure 
[?] shows that the temperature distribution of material, both dissipated and shocked, is accurately 
following self-similar evolution. The high resolution and standard resolution simulations of the 
n = 0, tc = 1 model also correspond well, suggesting that this distribution is insensitive to 
numerical resolution effects. As one would expect, the models with stronger cooling show an 
increased fraction of the baryons at cooler temperatures. As we increase the strength of the 
cooling in the n = models the distribution of gas temperatures becomes wider, developing a 
more extended low temperature tail that contains the bulk of the gas. In the n = +1 models the 
temperature distribution maintains a bell shape and simply shifts toward lower mean temperature 
with increased cooling. For tc < 200, the temperature distributions are also less sensitive to the 
cooling amplitude in the n = +1 models than in the n = models. 

For both n = and n = +1, we find that some gas substantially hotter than Th e is 
incorporated into the I = 0.2 Ax p groups, especially in the models without cooling. While this hot 
gas must be shocked, it may represent a different population (presumably hot, diffuse gas on the 
outskirts of these structures) from the interior, cooler gas at T <J Th e . This extra population of hot 
gas may help explain the surprising behavior noted in our discussion of Figure ||, where we find 
that the emission weighted temperatures of objects identified with the smaller linking parameter 
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Fig. 7. — Average gas temperature distributions for groups in the n = (left column) and n = +1 
(right column) simulations. Groups are identified by applying the friends-of- friends algorithm with 
£ = 0.2 Ax p to the full particle distributions (baryons and dark matter). Temperatures are scaled 
to The = /J.m p (2kB)' 1 G0, where 9 is the total linear mass density of the group. The solid curves 
are measured for groups at a/aj = 1 and the dotted curves for groups ata/a/ = 0.4. In the n = 0, 
tc = 1 panel we also plot the high resolution, N = 2 x 256 2 simulation at a/aj = 1 with the 
short dashed lines and at a/aj = 0.4 with long dashes. Since all quantities in this plot are defined 
in a scale- free manner, the profiles measured at different times should lie on top of one another 
if these structures are correctly obeying the expected self-similar behavior. In order to calculate 
these average distributions, the baryon mass in each group is binned as a function of T/Th c , and 
the individual distributions for groups within a given range of (0.5 9 n \ < 9 < 2.5 9 n \ for n = 0, 
and 0.2 9 n \ < 9 < 1.0 9 n \ for n = +1) are then averaged together. 
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actually scale better than those found with the large £, since the objects selected with the smaller 
linking parameter do not include this extra gas. 

In all cases the temperature distributions in Figure ^ are fairly continuous, while in typical, 
realistic cold dark matter (CDM) simulations (see, e.g., Katz, Hernquist, & Weinberg 1992 or 
Evrard, Summers, & Davis 1994) the temperature distribution of gas in such objects is more nearly 
bimodal, with dissipated gas having cooled efficiently to T ~ 10 4 K from shocked temperatures 
of T ~ 10 5 - 10 6 K. Though we do not show the results here, we have repeated our 2-D models 
using a realistic cooling law such as that used in Katz et al. (1992), and in such cases we recover 
the typical bimodal temperature distribution, with cooled, dense gas knots embedded in hot, low 
density gas halos. This test shows that it is our featureless power-law cooling relations that cause 
the temperature distributions to have the continuous form seen in Figure |7], not the difference in 
dimensionality. 

We are especially interested in the fraction of the baryons that cool to low temperatures 
following collapse and shock heating in dark matter halos. In full-scale simulations of galaxy 
formation, it is these radiatively cooled baryons that are converted into stars, creating the 
luminous regions of observable galaxies. If cosmological simulations are to correctly predict galaxy 
luminosity functions, star formation rates, and so forth, they must first correctly calculate the 
masses of the cold baryon concentrations that form the galaxies. 

We will turn to the self-similar scaling of the cooled baryon fractions shortly, but we first 
examine the dependence of this fraction on total mass at fixed time. Figure |8| plots {6c /6b) as a 
function of total group mass 9 at the final output time, where 9c /6b is the fraction of baryons in 
the group that have cooled to T < 0.005Th e , and the averages are computed in logarithmic bins 
of 9. We choose a low temperature threshold so that we probe only regions that have undergone 
radiative cooling, and as a result only the strong cooling models (tc = 1 and tc = 0.1) show up on 
the plot, along with two points for tc = 200, n = +1. The tc = 0.1 points (circles) show a clear 
trend of increasing 9c/9b with increasing 9, with a steeper correlation for n = +1 than for n = 0. 
The tc = 1 points (triangles) have larger scatter but appear to show a similar trend, especially 
with the higher dynamic range afforded by the 2 x 256 2 , n = simulation. 

These trends are expected given the temperature dependence of our cooling functions. The 
friends-of- friends algorithm selects groups with similar overdensity at all 9, and equation (|2|) 
predicts that the characteristic temperature of a group in the absence of cooling is proportional 
to 9. The cooling time is proportional to pT/A(T) oc T 1_/3 for groups of similar p, where (5 is the 
index of the cooling function given by equation (|l5| ) . For (3 > 1 this argument implies a shorter 
cooling time in more massive groups. The lines in Figure || show the relations 



^ocf^ocT^oc^- 1 , 




i.e., a cooled baryon fraction inversely proportional to the cooling time. While the relation 
between cooled fraction and cooling time need not be this simple, these lines do seem to capture 
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the numerical results fairly well, explaining both the steep trend of cooling with mass and the 
marked difference in slope between n = and n = +1. 

Since a given group has gas at a range of densities and temperatures, even in the absence of 
cooling, our argument above makes the implicit assumption that groups of different mass have 
similar density and temperature profiles in scaled units before the onset of cooling: formally 
speaking, p/p and T/T^ e have the same functional dependence on r/r v ; r , where r v i r is the group 
virial radius. This type of similarity is not guaranteed by our choice of scale-free initial conditions, 
which only ensures self-similar time evolution. (For example, scale-free physics implies that a 
group at time t\ is similar to a 0* group at time t2, but it does not imply that a 0.50* group 
is similar to a 9* group at fixed time.) This assumption is justifiable in light of recent N-body 
experiments, such as those of Navarro, Frenk & White (1995, 1996, 1997) or Cole &; Lacey (1996), 
who find that collapsed structures in high-resolution N-body experiments tend to a universal 
density profile over roughly two orders of magnitude in mass. We have directly examined the 
objects in our experiments and verified that the density and temperature profiles are similar over 
the mass ranges we resolve. 
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Fig. 8. — Average fraction of baryon mass per group with T/T^ < 0.005 as a function of total 
group 9 at a = aj for the n = (left) and n = +1 (right) simulations. The lines show the expected 
relations 9c/9b oc # 3//2 (n = 0) and 9q/9b oc 9 9 / 2 (n = +1), based on the assumptions that objects 
of differing masses are similar and that the cooled mass is inversely proportional to cooling time, 
as described in the text. Point types as defined in Figure ^: stellated triangles for models with 
cooling time tc = 200, open triangles tc = 1, circles tc = 0.1, and filled triangles the N = 2 x 256 2 , 
tc = 1 model. 
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Fig. 9. — Average fraction of baryon mass in groups with T/T^ e < x as a function of expansion 
factor for the (a) n = and (b) n = +1 models. At each expansion factor, the fraction of baryon 
mass 9c /Ob with T/T^ e < x is calculated for each group (where x is given in the right side row 
labels), then averaged for all groups in a given range of 6 (given in column headings). Since all 
quantities in this figure are defined in a scale-free manner, the cooled baryon mass fraction is 
expected to remain fixed throughout the evolution, as shown by the lines. Point and line types as 
defined in Figure |3|: stellated triangles/dotted lines for models with cooling time tc = 200, open 
triangles/short dashes tc = 1, circles/long dashes tc = 0.1, and filled triangles/dot-dash lines the 
jV = 2 x 256 2 , tc = 1 model. 
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Figure ^ directly addresses the central issue of our investigation, the scaling of the cooled 
baryon mass fraction. We take each object identified at a given output time and calculate the 
fraction of the baryon mass that has cooled below a fixed fraction of the hydrostatic equilibrium 
temperature, T/T^ e < x- We then compute the average of this mass fraction for all groups in a 
given range of group mass and plot these values as a function of the expansion factor. 

The upper left panel of Figure |9|a shows the average fraction of gas with T/T\ 1C < 0.2, for 
groups with mass in the range 0.16> n i to 6 n \. As expected, this fraction is higher for models 
with shorter cooling times: 15%, 75%, and 90% for tc = 200, 1, and 0.1, respectively. While 
self-similarity does not tell us what this fraction should be for a given model, it does tell us 
that the fraction should be independent of time, since the ratio T/Th e is dimensionless and the 
mass range is scaled in terms of 9 a \. The simulation results agree with this analytic scaling law 
remarkably well, with cooled gas fractions that are independent of expansion factor in the range 
log a/df € [—0.4,0]. The standard and high resolution simulations of the tc = 1 model are in 
in nearly perfect agreement for log a/aj > —0.4, and the high resolution simulation extends the 
scaling back to log a/cif = —0.8. If numerical resolution effects were important, we would expect 
to find a difference between these two simulations, and we would expect the cooled gas fraction 
to increase with time in any given simulation as groups in the range 0.19 n \ — 6 n \ are resolved by 
progressively more and more particles. 

The lower left panel of Figure ^a shows the average fraction of mass in these groups with 
T/Th e < 0.05. While less gas has cooled below this lower temperature threshold, the cooled 
gas fractions still show the analytically predicted scaling. The more massive groups have 
more cooled gas (right column of Figure ||a), as expected on the basis of Figure || and the 
accompanying discussion. The cooled gas fractions scale well over a larger range of expansion 
factor, log a/af G [—0.8,0], because these more massive groups are well resolved at earlier times. 

Figure shows the corresponding results for the n = +1 models. We select different mass 
ranges that correspond better to the numerical limits of the n = +1 simulations (see Figure |l|), 
and we adopt lower temperature thresholds because cooling is stronger in the n = +1 models (see 
Figure ^). The cooled baryon fractions of the n = +1 groups scale well for log a/aj G [—0.6,0] in 
the mass range O.O40 n i - 0.26» nl and for loga/a/ G [-1-0,0] for 0.20^ - 1.0<9 nl . 

In the simulations of Paper I, with no cooling, we found poor scaling of group Bremsstrahlung 
luminosities. We attributed this failure to poor scaling of the group's core densities — as groups 
are resolved by progressively more particles, the estimated central densities increase because of 
better resolution. Sensitivity of the gas density to numerical resolution has also been noted in a 
number of previous investigations, including Kang et al. (1994), Anninos & Norman (1996), Frenk 
et al. (1996), Weinberg, Hernquist, &: Katz (1997), and Owen &: Villumsen (1997), to name merely 
a few. 

Since radiative cooling is itself sensitive to gas density (A tx p 2 ), we initially expected that 
it would be difficult to find self-similar scaling in simulations with radiative cooling. How, then, 
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do we account for the nearly perfect scaling of the cooled baryon fractions in Figure |9|? The key 
to the explanation is that the cooling instability is effectively a threshold phenomenon. If the 
numerical resolution in an object is insufficient, so that the cooling time at the highest resolved 
density is longer than the dynamical time, then radiative cooling is almost entirely suppressed. 
Once the resolved density passes a critical threshold, however, thermal instability sets in, and 
the gas rapidly cools, loses pressure support, and converges toward the physical solution. The 
fact that these simulations do not form significant amounts of radiatively dissipated gas until the 
larger objects begin to demonstrate the expected self-similar scalings (as noted in our discussion of 
the mass and temperature scaling tests in Figures ||] and f|), indicating that the physical solution 
is beginning to be followed correctly, supports this interpretation. The agreement between the 
standard and high resolution simulations of the n = 0, tc = 1 model and the successful scaling 
of the mass function of dissipated objects (Figure ||), the temperature distributions in collapsed 
groups (Figure ||) , and the cooled baryon fractions (Figure ^) all indicate that the simulations have 
converged to a physical solution that is, at least for these measures, insensitive to their numerical 
resolution. 

7. Summary and Conclusions 

We analyze a set of 2-D hydrodynamic simulations designed to study the self-similar evolution 
of hierarchical structure in the presence of radiative cooling. The radiative cooling law for a 
primordial H/He plasma would introduce a preferred timescale in these models, so we instead 
use artificial cooling laws that maintain the scale-free nature of the physics: A(T) oc T@ where 
P is a function (eq. (T^j) of the spectral index n of the initial density fluctuations. We simulate 
eight distinct physical models, with n = and n = +1 initial power spectra and four different 
amplitudes of the corresponding radiative cooling laws, ranging from no cooling to very strong 
cooling. For each of these models we evolve a simulation with 128 2 gas and 128 2 dark matter 
particles, using ASPH to model the hydrodynamics and a Particle-Mesh technique to follow the 
gravitational interactions. We also repeat one model with 256 2 gas and dark matter particles, to 
directly assess any numerical resolution dependencies. 

Because both the physics and initial conditions of these models are scale-free, their physical 
properties must evolve self-similarly in time, in accord with analytic scaling laws. These rigorous 
analytic predictions provide a stringent test of our numerical methods, since numerical parameters 
like particle mass and force resolution stay fixed as the system evolves and will therefore act 
to break self-similar behavior if they limit the physical accuracy of the calculation. We test 
for self-similar evolution by identifying objects as groups of particles within a given overdensity 
contour at various times, measuring global properties characteristic of these structures, and 
examining the evolution of these quantities over a range of output times. In general we find 
excellent agreement between the analytically predicted and numerically measured scalings for 
the masses, temperatures, Bremsstrahlung luminosities, and (to a lesser extent) gas densities 
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over the range of expansion factors that we would naively expect to be accessible based on the 
simple resolution arguments presented in §|4|. Our most significant result, demonstrated directly in 
Figure |9| and elaborated in Figures |2| and ||, is that the fraction of baryonic material that cools in 
collapsed objects of specified mass follows the expected scaling with remarkable accuracy 

One of the impressive successes of hydrodynamic cosmological simulations with CDM initial 
conditions and realistic cooling laws is that they produce dense clumps of cold gas with masses, 
sizes, and overdensities comparable to the luminous regions of observed galaxies (Katz et al. 1992, 
1996; Evrard et al. 1994; Summers, Davis, & Evrard 1995; Frenk et al. 1996; Weinberg et al. 1997; 
Pearce 1998). With reasonable prescriptions for galactic scale star formation, these objects are 
converted into dense, tightly bound clumps of stars and cold gas. The resulting stellar masses are 
not sensitive to the parameters of these prescriptions, at least within some range, because the star 
formation rate is governed mainly by the rate at which gas cools and condenses onto the central 
object (Katz et al. 1996; Pearce 1998). If numerical simulations are to provide accurate predictions 
of quantities like the galaxy luminosity function, star formation rates in high redshift systems, 
or even the bias between galaxies and mass, then they must accurately follow the gravitational 
collapse, shock heating, and subsequent dissipation and condensation of baryons into these dense 
systems. 

Our results provide encouraging evidence that cosmological simulation methods can indeed 
rise to this challenge. Specifically, they show that 2-D calculations that resolve individual systems 
with a few dozen to a few hundred particles correctly follow the formation of dissipated objects. 
They also suggest that computing the galaxy baryon mass function may be an easier problem 
than computing the cluster X-ray luminosity function despite the higher densities of the objects 
in question and the greater complexity of the physics (with radiative cooling playing a larger 
role). The p 2 dependence of Bremsstrahlung emissivity implies that a simulation must resolve an 
object's central density in order to compute its X-ray luminosity accurately. However, thermal 
instability is a threshold phenomenon, and once a simulation resolves an object's density at the 
cooling radius (where the cooling time equals the dynamical time), it can compute the dissipated 
baryon mass with reasonable accuracy. While it is clearly important to repeat the self-similar 
evolution tests with 3-D calculations, we see no reason that the behavior in three dimensions 
should be fundamentally different. 

For nearly two decades, cosmological N-body simulations have provided an indispensable 
tool for predicting the large scale distribution of dark matter in different cosmological models. It 
appears that the more ambitious goal of predicting the properties and distribution of galaxies with 
hydrodynamic simulations is now well within reach. 

This work was supported by NASA grants NAG5-2882 and NAG5-3111. JMO's work was 
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